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Electrodynamics becomes nonlinear and permits the self-interaction of fields when the quantised 
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I. INTRODUCTION 

Soon after the formulation of relativistic quantum 
mechanics, it became clear that the propagation of 
light through the vacuum would be modified due to the 
polarisability of virtual electron-positron pairs [1-4]. 
Heisenberg and Euler derived the Lagrangian of an 
effective description of this interaction for constant 
fields [5], which was later rederived by Schwinger [6]. 
Derivative expansions of this effective interaction [7-9] 
and numerical worldline calculations [10] imply that 
“constant” is to be taken with respect to the Compton 
time h/m(? for electron mass m. This suggests a good 
approximation of the effect for time-dependent fields 
with a much longer period than the Compton time is 
to simply insert them in place of the constant fields 
in the Heisenberg-Euler Lagrangian. In particular, 
the polarised vacuum supports the phenomenon of 
self-interaction when two electromagnetic waves couple 
via virtual electron-positron pairs and the principle of 
superposition no longer holds. 

There have been several studies of the consequences 
of this self-interaction. Lutzky and Toll [11] showed 
that if the field invariant Q = —FF*/4E^j. = E • B = 0 
where F, F* are the Faraday tensor and its dual, 
Ecr = m^c^/he = 1.3 • 10^®Vcm“^ is the so-called 
“critical” field, e > 0 is the charge of a positron and 
E and B the total electric and magnetic fields in units 
of the critical field, a current that depends nonlinearly 
on the invariant F = —F'^/AE^^ = {E^ — B^)/2 leads 
to the generation of an electromagnetic discontinuity or 
“shock”. After identifying an application in magnetised 
neutron stars, shocks were analysed in a constant 
magnetic field background using a first- [12], second- [13] 


* patrick.boehl@physik. uni-muenchen. de 
t hartmut.ruhl@physik.uni-muenchen.de 
t b.king@plymouth.ac.uk 


and several- [14, 15] order weak-field expansion of the 
Heisenberg-Euler Lagrangian with an all-order analysis 
performed by Bialynicka-Birula [16]. An astrophysical 
environment was further modelled by introducing 
nonlinear vacuum effects into equations of relativistic 
magnetohydrodynamics [15] and into a dusty plasma [17]. 

In the current article we analyse a pump-probe set-up 
of having a linearly-polarised oscillating plane wave 
(probe) counterpropagate through a linearly-polarised 
constant crossed and stronger plane wave background 
(pump). Of particular interest will be the two cases 
of having parallel or perpendicular probe and pump 
polarisations. Observables are expressed in terms of 
the electric and magnetic fields to aid comparison with 
numerical simulation. 

Unlike in classical electrodynamics where a super¬ 
position of solutions to the wave equation is also a 
solution, when the existence of charged virtual electron- 
positron vacuum states is included, the principle of 
superposition is no longer valid. A consequence of using 
the Heisenberg-Euler Lagrangian is that a non-trivial 
vacuum “current” appears in Maxwell’s equations, which 
disappears in the classical limit h ^ 0. If the electro¬ 
magnetic fields are not very weak E ^ y/a{tiuj/mc^)'^ 
[18] where a ~ 1/137 is the fine-structure constant and 
the field frequency is w (corresponding to an intensity 
much greater than 10® Wcm“^ for an optical laser), they 
can be regarded as classical. When this is the case, the 
methods of classical electrodynamics can be used to 
solve Maxwell’s equations with the vacuum current. 

For fields much weaker than critical, the interaction 
with the virtual electron-positron pairs of the vacuum 
permits 2n-wave mixing for integer n > 1 such as four- 
and six-wave mixing, as demonstrated in Fig. 1. One 
can make an analogy with nonlinear optics, in which the 
polarisation P of an optical material can depend upon 
higher powers of the electric field [19], which are de¬ 
scribed using different orders of the susceptibility tensor 
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FIG. 1. An illustration of the weak-field expansion of the 
vacuum polarisation diagram. 


^0). 

P^ = X^]^E,+x^,lE,Ek + ..., ( 1 ) 

and analogously for the magnetisation M. Being 
a relativistic effect however, the magnetic and the 
electric field appear in the vacuum polarisation and 
magnetisation on an equal footing. For weak fields 
and propagation lengths shorter than the scattering 
length, four-wave mixing is the most probable vacuum 
polarisation process for colliding plane waves (with the 
exception of certain special field geometries). This is 
often compared to the optical Kerr effect [20], but in a 
pump-probe experiment in which the probe oscillates 
much quicker than the pump field, the steepening of 
the carrier wave [21] and not the envelope [22] is stronger. 

If the fields’ spacetime extent is much larger than a 
single scattering length, multiple 2n-wave mixing can 
occur, in which the change in field due to wave mixing 
influences further changes due to wave mixing, with 
each mixing event involving a potentially different n. 
Usually it is assumed that the probability for multiple 
mixing events is much lower than single mixing events, 
and multiple events are neglected. However, if the 
extent of the field is large enough, this hierarchy can 
be broken and it can become more probable that 
multiple mixing events occur than a single mixing 
event so that all orders of wave mixing events have 
to be taken into account. With the generation of a 
large number of higher harmonics, the shape of the 
electromagnetic plane waves will also change and this 
leads to the possibility of shock wave generation. In 
the “shock regime”, as all orders of wave mixing can 
play a role in the generation of the spectrum, the 
spectrum is expected to be qualitatively different to the 
perturbative case of having only a single mixing event, 
where four-wave mixing is the most probable and higher 
harmonics are exponentially suppressed. Such a type of 
shock generation is also known from nonlinear optics [23]. 

In contrast to this, the weak-field expansion in Fig. 
1 suggests that high harmonic generation can also 
occur through single scattering events that involve large 
numbers of photons. The likelihood of this happening 
increases with the field strength of slowly-varying weak 
fields. This type of vacuum high harmonic generation 
has been investigated using the full polarisation operator 
in [24, 25] and using the lowest order of the weak-field 
expansion in [26-29]. A highlight of the current article 


is the first investigation of vacuum high harmonic 
generation in plane wave fields in what we call the shock 
regime, where the probe propagation length is much 
larger than the mean scattering length. It will be shown 
that in certain parameter regimes, this can be a much 
more efficient high harmonic generation mechanism. 

The aims of this work are: i) to investigate vacuum 
high harmonic generation in the collision of plane waves 
that are weaker than critical, for the case that the 
fields’ spacetime extent is much larger than the mean 
scattering length; ii) to show that the higher harmonics 
are accompanied by an electromagnetic shock due to the 
polarised vacuum; iii) to investigate the dependency of 
this shock on the colliding fields’ mutual linear polarisa¬ 
tion; iv) to comment on the similarities and differences 
of high harmonic generation in laser-irradiated plasmas. 

We begin with a derivation of the modified Maxwell 
and wave equations (Sec. II), summarise the analyti¬ 
cal method (Sec. Ill) and the numerical method used 
in computational simulation (Sec. IV) before analysing 
higher harmonic generation with just four-wave mixing 
(Sec. V), just six-wave mixing (Sec. VI) and both four- 
and six-wave mixing (Sec. VII). We then discuss the 
results, compare with high harmonic generation from os¬ 
cillating plasmas (Sec. VIII) and conclude (Sec. IX). 


II. MODIFIED ELECTROMAGNETIC WAVE 
PROPAGATION 


The Heisenberg-Euler Lagrangian can be written [30] 

[ 6 ] 


>Che = 


Stt^ 



s^ab cot as coth 6s — I 



( 2 ) 


(we have set here and throughout h = c = 1 unless they 
explicitly occur), where the secular invariants a and 6 are 
given by: 


a = 




1/2 


6 = 




-, 1/2 


and we recall that electric and magnetic fields are in units 
of the critical field E^r- Applying the Euler-Lagrange 
equations to £ = £mw + £he, where £mw = nP{E‘^ — 
B'^)/8Tra is the classical Maxwell Lagrangian, gives the 
modified Maxwell equations: 


= 0 (3) 

(1 -b Cl) + C2 + C3 F*>^''d^{FF*) 


Ca 


F*t^^df,F^ + Ff^'^d^iFF*) =0 (4) 


and the general expressions for the coefficients Ci are 
given in App. A. Expressing these equations in electric 
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and magnetic fields, we acquire: 

V A E + 5tB = 0 (5) 

VAB-9tE = J[E,B] (6) 


J[E, B] = [Cl (9tE - V A B) + (CzE + C 4 B) dtF^ 

+ (C 2 B - C 4 E) A 

+ (C 4 B - C 3 E) A V(FE*) ^ ’ 

+ (C3B + C4E)5t(FE*)]. 

The current J in Maxwell’s equations is related to the 
corresponding source T in the wave equation for the 
electric field via J = T. 


We restrict our analysis to the case when E <C 1, for 
two reasons. First, it allows us to neglect the creation of 
real electron-positron pairs, as the probability of vacuum 
pair production in a volume equal to the reduced 
Compton wavelength X = h/mc cubed in the Compton 
time A/c is P = exp{—n/E)/A tt^ [6], which is heavily 
suppressed for E <C 1. Second, it permits a perturbative 
expansion in E, the so-called “weak-field expansion”, of 
the Heisenberg-Euler Lagrangian. 


Although all electromagnetic fields are classical, it is 
useful to envisage the corresponding quantum process in¬ 
volving photons and this is depicted for the weak-field 
expansion of the vacuum polarisation operator in Fig. 1. 
(Indeed, it has been shown that the leading-order term of 
the weak-field expansion agrees with the direct calcula¬ 
tion of the four-photon box diagram in the low-frequency 
limit fiLo <C mc^ [31]-) The weak-held expansion of Eq. 
(2) for E 1 is then: 


Ehe 

Cl 

C 2 

C 3 


a 


7(E- B)^ 


/A 
47r 

^(E^-E^) 2(e2-E2)%13(E-B) 


( 8 ) 

(9) 

? 

■( 10 ) 



3 (E^ - E^)'^ + 22 (E^ 
19 (E • B)^l , 


e2)^(E-B)^ 

( 11 ) 


where \ii = “/qott, ^2 = /ts = 4a/9457r (although 

a occurs in the denominator in Eq. (8), as helds are in 
units of the critical held, when H ^ 0, Ehe —^ 0). The 
coefhcients Ci in Eq. (4) that follow from £1 and £2 are 
given by Eqs. (A5-A10) in App. A. 

In the scenario we consider, the initial electric held is 
E(°)((/3p, ips) = Ep°^((/jp) + e(°^((/?s) and the initial probe 


and strong electric waves are given by: 

E^°H£p) = cosv3p (12) 

Ef(:^,) =e,£,Rect (^1^^ , (13) 

where the rectangular function Rect((^/4>) = 

9{(p + 4>/2) — — <i>/2) and d(-) the Heaviside 

function [32], ipp = kpX = ujpX~, (ps = kgX = ujgX'^, 
x^ = t ± z, ^p = u!pTp, = uJsTg with the probe 
and strong held polarisation vectors e^, Eg obeying 
— 1, — 1; kp ■ Sp — 0, kg ■ £g — 0 and 

the probe pulse is assumed to be much weaker than 
the strong background £p ^ Eg. Initially, the probe 
and strong helds are well separated: limt_>._oo E,G = 0. 
We dehne the orthonormal polarisation vectors (ell,e''‘) 
where = £p dehnes “parallel” polarisation, and £"*■ 
“perpendicular” polarisation with -ell = 0, -kp = 0. 

Since the vacuum current is a function of the rela¬ 
tivistic invariants T = (E^ — E^)/2 and C/ = E • B, for 
single plane waves, there is no effect on propagation due 
to vacuum polarisation [6, 33, 34]. Therefore, the only 
contributions will come from cross-terms between the 
probe and strong held. As the weak-held expansion is 
an expansion in powers of E and G, for our scenario, 
each order scales as £„ ~ (EsEp)"+^. 

The initial probe Ep°^ and strong eI''^ helds satisfy the 
classical vacuum wave equation independently: 

□ e[,°) = 0 , □ Ei°) = 0, 

where □ = — V^. The effect of the polarised vac¬ 

uum can be included with a source term T = T[E, B] oc¬ 
curring on the right-hand side of the wave equation. We 
will assume that solutions to this equation are also plane 
waves propagating along the same axis as the pump and 
probe waves. This allows us to write T = T[E]. Since a 
single plane wave cannot polarise the vacuum [33, 34]: 

T[Ep] = 0 , T[E,] = 0. 

However, since two counterpropagating plane waves can 
polarise the vacuum, the wave equation we will solve is: 

□ (Ep + E,) = T[Ep + E,]. (14) 

In particular, we are interested in solutions which in¬ 
clude the self-action of the probe that lead to a plasma¬ 
like vacuum instability and corresponding electromag¬ 
netic shock. Eq. (14) will be solved in two ways. Eirst, 
the scattered probe will be solved for using an analyti¬ 
cal method based on an iterative procedure that ignores 
changes to the stronger background: 

□ Efo-Hi) ^ e(o)]^ (^5) 

Second, Eq. (14) will be solved consistently in a 
numerical simulation that uses tools based on the 
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pseudocharacteristic method of lines, which are applied 
to the corresponding Maxwell equations. In this way, 
the “asymptotic” state of the probe field after it has 
passed through the strong Held and T r; 0 (in contrast to 
the “overlap” dynamics when T 7 ^ 0 [28]) will be studied. 


analysis. An argument for neglecting eight-photon scat¬ 
tering will be forthcoming. 

III. ANALYTICAL METHOD 


As we are considering the collision of counter- 
propagating plane waves, the general Maxwell’s equa¬ 
tions in Eqs. (5) and ( 6 ) reduce to one spatial (z) and one 
temporal (t) dimension. To determine which terms in the 
full weak-field expansion for the current Eq. ( 8 ) should 
be considered when calculating high harmonic genera¬ 
tion, we employ the following scaling argument. As ex¬ 
plained in [28], the change in the field due to interaction 
with the vacuum that propagates with the probe (“for¬ 
ward” scattering) is 


To solve the inhomogeneous wave equation 

[a,^-a2]Ep = T[Ep + E(°)], (18) 

we employ an iterative ansatz: 

E^”+i) =EW-f AE^, (19) 

where 


AEp(x-)=/ — J(t' = x- + z',z'), (16) 

J —00 ^ 


where the vacuum current is: 


J — ^ Jj ; Ji _ 47r 


i=l 


k„ A dzM^ + StPi 


(17) 


kp = kp/jkpj and the dimensionless vacuum polarisation 
Pi = 9£i/9E and magnetisation = 9£i/9B (as used 
in e.g. [16] or [35]). The forward-scattered signal is zero 
if the vectorial part of Pi or Mi is from the probe field. 
As already explained, £„ ^ (Es£'p)"+^, but in the wave 
equation that results from this, the vacuum current con¬ 
tains derivatives with respect to Eg and Ep. Since the 
current containing the derivative with respect to Eg van¬ 
ishes for forward scattering in plane waves [28, 36], we 
see that the remaining current and hence the scattered 
field Jn oc . The integration over z' is over the 

strong field and so contributes a factor Tg and the dif¬ 
ferentials in Eq. (17) contribute approximately a factor 
Wp, so that one can estimate AEp^^ oc $, for 

4) = LOpTg. Since we assume E <g; 1 , and since we are 
interested in the case when the change in the probe is 
of the same order as the probe field and self-interaction 
becomes important, we require $ ^ 1. We also note 
that the coefficients /r„ diverge with n because the weak- 
field expansion is asymptotic (see e.g. [37]), so we do 
not expect the series can be truncated for arbitrarily 
large n and still yield a useful approximation. Although 
purely four-photon scattering does allow the generation 
of higher harmonics in this set-up, this first occurs for 

double four-photon scattering. The contribution from 

i2') 

this twice-iterated process appears in AEp ' and scales as 
oc (/ii)^£’^£’p4), which when compared to the leading con¬ 
tribution to second harmonic generation from six-photon 
scattering in AEp^'^ oc ^ 2 ^ 5 ^p^h, is suppressed by a factor 
(iJiiY/IJi 2 “C 1. Therefore, when considering higher har¬ 
monic generation along the probe propagation axis in the 
regime E 1 , 4’ ^ 1 , the leading contribution originates 
from six-photon scattering. In Sec. V, this simple scal¬ 
ing argument will be seen to agree with the full numerical 


and in general 

00 

T(")(t,z) = [e("-i)(^p) +E(0)((^,)' 


i=l 


where the subscript i is the order of the weak-field ex¬ 
pansion and the retarded Green’s function is [38]: 

G{t.z) = ^o{t)e - [z]^ , 


for refractive index n. If n = 1, one acquires Eq. (16), 
where = T^^'>{t,z). These equations can be 

iterated to calculate the generation of higher harmonics 
due to multiple scattering as outlined in the introduction. 
Within this analytical approach, we assume uipTp ^ I 
and uipTg ^ 1 , so that the derivative of the probe and 
background envelopes can be neglected with respect to 
the derivative of the oscillating part of the probe in J. 
When studying the generation of higher harmonics, we 
will be particularly interested in taking 


T(-)(t,z) = T2 [e("-i)(^p) -f Ei°)(^,) 


which corresponds to considering purely six-photon scat¬ 
tering (this will be further justified shortly). 

A diagrammatic approach is useful to understand the 
physical processes described by different iterations of the 
probe field Ep”\ Eirst, since we are interested in har¬ 
monic generation and since the background is constant, 
we suppress strong-field photon legs. Eurthermore, as 
the Heisenberg-Euler Lagrangian is “effective” in that all 
fermion dynamics have been integrated out, all vacuum 
loops are reduced to effective vertices. Then the diagram 
representing six-photon scattering, which is the leading 
order harmonic-generating process, is given in Eig. 2. 
The iterative ansatz in Eq. (19) is illustrated in Tab. I. 


The diagrammatic equation in Tab. I in some ways 
resembles the Schwinger-Dyson equation [39] but in this 
case the left-hand side is the self-consistent solution of 
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rtojp \ 

\ 0, ±.2ujp 

/ 

iojp / 


FIG. 2. In the left-hand diagram, ujj € 3tJs, 2a;p ± 
LjJa,2u}p ± 3cJs}. If the strong field is approximated as con¬ 
stant and the three strong-held photon legs are suppressed, 
in an effective approach, six-photon scattering of the probe 
can be represented as a triple interaction. The ± refer to 
incoming and outgoing photons respectively. 


E 


("+!) 

P 


e: 


(0) 


AE^”^ 



pi(o) . 

thp 

thp 


E 


( 2 ) . 
P 



TABLE I. Diagrammatic representation of the hrst iterations 
of the probe wave equation. 


the probe field at a particular order of iteration, and 
the double line on the right-hand side is where the 
scattered probe field from the previous order is applied. 
In Tab. I, it is shown how the number of diagrams 
rapidly increases with iteration order (as the square of 
the number in the previous order plus one, although 
many are equivalent). It also demonstrates that terms 
of a much higher perturbative order (number of vertices) 
are generated at a given iterative order (E*^") contains 
terms from the (2" — I)th perturbative order, but is only 
accurate to the nth perturbation order). 

On all the diagrams with at least one vertex, one 
leg is the scattered field and the rest are incoming or 
outgoing probe photons. An example is given in Fig. 


3 where the ± sign refers to the energy added to the 
system by incoming/outgoing photons. By summing 

i idp 

i cjp-q 

p 0, 2 ujp, 4 ujp 

i ujp -q 

i ujp 

FIG. 3. An example of the harmonics generated in the 
probe due to effective self-interaction in a slowly-varying back¬ 
ground. 


the series that occurs in lim„_>oo Ep"^, we will arrive at 
an analytical expression for the asymptotic probe field 
and in doing so identify a shock parameter that signifies 
when self-action effects become important. 


For the example of parallel probe and strong field po¬ 
larisation, the second iteration shown in Tab. I is; 


I]^^'>=epSpe Up 


“ (i) COS(pp 

“2^^^ - 3 g^^^\ips)cos‘iipp 

+2(^) g(^)(v2s)sin4(^p , (20) 


where v = 1^2 exp(—((pp/$p)^) and the shock parameter 
V 2 = 192^2^3^*^- functions of ifg describe how 

the particular term is generated during the passage of 
the probe through the strong background (all fields are 
classical) and originate from repeated integration of the 


interaction over co-ordinate. Here: 



/ dy Rect(y) 

/ — 00 


II 



/ dy Rect{y) g^^\y) 

II 



/ dy Rect(y) 

' —00 

9^^Hy) 


and these are plotted in Fig. 4. 

As mentioned in the introduction, we are mainly inter¬ 
ested in the asymptotic state of the probe: 

E(2)(^p)= hm Ef((^p,;^,), (22) 

where we note lim,p^_>.oo = 1, 

lim,p,,^oo5*'“nv?s) = V 2 and g^'^Ky^’s) = V3- 

As previously remarked, using this method, Ep"^ con¬ 
tains powers of v from 0 to 2" — 1 but is only accurate 
to 0{v^). We also note that the nth iteration generates 
harmonics from I to 2". A power series in v multiplies 
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FIG. 4. (Color online) A plot of how the functions describing 
how the occurrence of higher harmonics varies with probe 
propagation length. 


each harmonic so we can write a given iteration as: 

=ep£pe~y^) ^ ^a‘^^\v,ipf,)sin2jipp 


t=i 


+4”^i(r;, (/Js)cos(2j - l)(/7p .(23) 


Of most interest is the asymptotic state of the full solu¬ 
tion: 


Ep((/ 3 p) = lim lim E(”)((/Jp,(/Js), 

(ps—^oo n—)-oo ^ 

and we find that for the parallel set-up: 

lim lim a['^\v,(ps) = cLj{v) = 2 (— 1 ) 0 ^) ^ 

c/?s—>-oo n—>-oo jy 

(24) 

where [jj = floor(j) and Ji{-) is the Zth-order Bessel 
function of the first kind [40]. We note that the all-order 
solution Eq. (24) for a plane probe propagating through 
a constant crossed, parallel-polarised background, re¬ 
sembles the Fubini solution [41] for the propagation 
of lossless finite-amplitude planar acoustic waves in 
nonlinear media [42]. 

The all-order solution can be derived from a probe- 
dependent refractive index: n = 1 -|- dn 2 with 1/2 = i 5 n 2 $ 
where: 


Sn 2 {ips,‘Pp) = 192fi2E^{<Ps)Ep{(pp). (25) 


To justify when it is a good approximation to only 
consider six-photon scattering, let us consider first eight- 
photon scattering. The shock-parameter for eight-photon 
scattering is 1^3 = 1536^3fgfp4>. In order that this is 
much less than j^ 2 , we require £s£p <C 2 / 32 , and since 
fs -C 1 and fp -C 1, this is fulfilled. Therefore the indi¬ 
vidual effect of the next higher-order terms in the weak- 
field expansion should be negligible. In contrast, the im¬ 
portance of four-photon scattering can be quantified by 
the parameter vi = $ but this corresponds to the 

process of one incoming and one outgoing photon from 
a scattering event and therefore will not contribute di¬ 
rectly to harmonic generation. Nevertheless, it does lead 
to a refractive index alteration, which in combination 
with multiple six-photon scattering, could potentially in¬ 
fluence the generated spectrum. To ignore this in our 
analysis would require £s£p S> which is not fulfilled. 
To explore this point, the simulation results are split into 
three cases: i) purely four-photon scattering, ii) purely 
six-photon scattering and iii) both four- and six- photon 
scattering. We consider two polarisation scenarios: the 
“parallel set-up” and the “perpendicular set-up”, which 
refer to the initial strong field polarisation being in the 
ell and mode respectively. In the parallel set-up, we 
will find all harmonics are generated in the parallel po¬ 
larisation mode ell, whereas for the perpendicular set-up, 
each odd harmonic will be generated in a perpendicular 
mode e^ and each even harmonic in a parallel one ell. 


IV. NUMERICAL METHOD 

For the scenario of two colliding plane wave pulses, the 
modified Maxwell equations in Eqs. (5) and ( 6 ) can be 
written in matrix form: 

(l 4 + X)atf+(Q + Y)a,f = 0, (27) 

where f = {E^, Ey, B^, By)"^, I 4 is the identity matrix 
in four dimensions, Q = adiag(l, —1, —1,1) is an anti¬ 
diagonal matrix and X and Y are the perturbations due 
to vacuum interaction given in a general form in App. 
B. 


Our numerical method, which was first employed 
by the authors in [28] and will be explained in more 
detail in the following, is based on inverting the matrix 
(I 4 -I- X) to convert Eq. (27) to a system of ordinary 
differential equations (ODEs), discretising in space using 
the “pseudocharacteristic method of lines” (PCMOL) 
[43] and integrating the equations of motion using the 
ODE solver CVODE [44]. 


So the scattered probe field due to just six-photon scat¬ 
tering can be written: 

Ep((^p) = E^°^ {ipp + V2[Ep{ipp)]). (26) 


Our analysis is valid when £s,£p ^ 1 and the sin¬ 
gle parameter relevant to high harmonic generation in 
the shock regime that depends on the field strength is 
V 2 = 192pL2£l£p^. We wish to simulate the occurrence 
of a shock wave, for which 1/2 1 , implying must 

be very large in order to compensate for the weak field 
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strengths. However, a large $ is computational expen¬ 
sive to simulate. To compare analytical and numerical 
results, we will therefore extrapolate the theoretical re¬ 
sult to values ot Eg ^ 1, allowing a simulation for smaller 
$ to be performed, with the condition that the physical 
prediction is only valid for a particular value of V 2 when 
fs <C 1. For this reason, we will often quote simulation 
parameters in terms of shock parameters rather than ab¬ 
solute field strengths and spatial extensions. 


A. Linear case 


Let us first consider (27) with X = Y = 0, which is 
the —>■ 0 limit. This system is hyperbolic [45], which 
means that we can find a basis u := S f such that the 
matrix A = SQS~^ = diag(—1, —1,1,1) is diagonal with 
real eigenvalues: 


S 


1 


-10 0 1 

0 110 
10 0 1 

0-110 


u: 



By - Ex 

Ey-\-Bx 

Ex-\-By 

Bx-Ey 


(28) 


In this new basis, we have an uncoupled system of ad- 
vection equations: 


of ODEs u'(t) = g[u(t),t], where g[u(t),t] = — ADu, 
with the 4x4 matrix A being mapped onto a AN x 47V 
dimensional block-diagonal one, A = Ijv 0 A (0 is the 
Kronecker product [47]) and D being the 47V x 47V 
matrix representing the biased differencing explained 
above. For the detailed action of D on u see Appendix C. 

The initial conditions are set up in f, the system is 
integrated in u using CVODE and transformed back for 
output. CVODE is an ODE-solver that offers variable- 
order, variable-step multi-step methods. Initially, we 
supply the “right-hand-side function” g[u(t),t] as above. 
Since both the linear and nonlinear cases are non-stiff 
(no rapidly-damped modes are expected), we apply the 
Adams-Moulton-Methods together with the variational 
method to solve the resulting linear system. This pro¬ 
vides higher accuracy with less computational effort com¬ 
pared to the offered Newton iterations, since neither ap¬ 
proximations nor an analytical expression for the Jaco¬ 
bian have to be provided. We always use the parallel im¬ 
plementation of CVODE together with “extended” (long 
double) precision. 


B. Nonlinear case 


dtu{t, z) -I- A dzu{t, z) = 0. 

The diagonal elements Xi of A are called the “charac¬ 
teristic speeds" the system, where = ±1 corresponds 
to a component travelling along the characteristics 
with the speed of light. We proceed by introduc¬ 
ing a co-located grid for the components Ui with TV 
grid points. The field components Ui on the grid are 
arranged blockwise in a large 47V-dimensional vector 
u = (... ...), where = Ui{lAz) and 

0 < I < N is the index of the grid point. The PCMOL 
uses biased differencing for each component Ui according 
to the sign of the corresponding characteristic speed Xi, 
where the component Ui with > 0 {Xi < 0) is thereby 
differentiated using backward (forward) finite differences 
using fourth-order accuracy. In [46] it is argued that 
this biased differencing using five-point-stencils is an 
effective fixed grid method for first order hyperbolic 
partial differential equations because it shows a good 
balance between introducing minimal numerical diffusion 
and oscillations in the solution where steep gradients 
are present. The derivatives at the boundary are also 
approximated using only field values inside the box. 
Instead of transforming the system back to f (the tilde 
in this section indicates the discretised version on the 
grid), which is normally done in the PCMOL, the system 
is solved for u. This has the advantage of having open 
boundary conditions since the components Ui are only 
allowed to flow in one direction. If we take the system 
to be of size L and a spatial resolution of TV grid points, 
then distance is measured in units of Az = ^/(N-i), 
where TV — I corresponds to the boundary conditions 
being taken into account. We are left with a system 


By discretising the full nonlinear system (27), the 
matrices X and Y also become 47V x 47V dimensional. 
The system then can also be brought into ODE form 
= g[u(7),7] by inverting the matrix (l 4 ]v + X). 
Since X depends only on the field components, the full 
matrix can be written as X = Ijv 0 X* and the upper 
index denotes the former 4x4 matrix X at grid point 
1. This can be used to reduce the inversion of X to TV 
times the inversion of a 4 x 4 matrix. The structure of 
X* allows us to rewrite X* as X* = G H* with 

/I 0\ 

^ ^ 1 ^12 ^13 ^14 

0 0 I ’ v 2:2 i 2:22 2:23 2:24 

Vo oj 


where the are the values of the non-vanishing matrix 
elements of X given in App. B, evaluated at position 1. 
Then we can apply the Woodbury Formula [48], 

(I4 + X')-i = I4 - G(l2 + H'G)^1 h' , 


to further reduce the inversion to one of the 2 x 2 matrix 


(12 + H'G) = 


'-11 


I 


‘-12 

-x \2 


This is performed for all grid points using an LU- 
factorisation at each evaluation of the function g[u(t), t]. 


For the parameters considered, the nonlinear correc¬ 
tions X and Y do not change the signs of the character¬ 
istic speeds, so we use the same biased differencing as in 



the linear case. The nonlinear ODE-systein is then given 
by 

u'(t) = -S(l4Ar + X)-i(Q + Y)S-^Du . 

where S = lAr(8)S, Q = lAr(8)Q and Y = l^r 0 Y* in 
analogy to X. All fields are normalised by Ecr- The pa¬ 
rameters for CVODE are the same as in the linear case. 
The signals are analysed under the assumption w = |k| 
using a spatial Fourier Transform in Wolfram Mathemat- 
ica [49]. 


C. Simulational Setup 

Recalling the form of the probe and strong pulses (Eqs. 
(12) and (13)), we consider a Gaussian probe pulse with 
base frequency ujp and a “constant” strong pulse. We 
consider the two cases of parallel and perpendicular set¬ 
ups which are characterised by parallel and perpendicu¬ 
lar polarisations respectively of the probe and the strong 
pulse with Ep • = 1 and Ep ■ Eg = 0. The rectangu¬ 

lar shape of the strong pulse is approximated using a 
mirrored Fermi-Dirac distribution in the simulation box. 
The function FD(j/) is given by 



FIG. 5. (Color online) The simulational set-up for the col¬ 
lision of a Gaussian probe pulse with a “constant” pulse of 
various strengths (indicated by different line styles) and iden¬ 
tical polarisation. The size of the system is taken to be 
3.2 • 10“^ cm. 


with the analytical model (see Eqs. (16) and (21)): 


FD(y) 


1 


1 -I- exp( 




(29) 


The parameters Zf, and Zm play the role of the “tem¬ 
perature” and “chemical potential”, controlling the 
steepness and width of the strong pulse. Typical values 
are Zf, = 5 • 10“®cm and Zm = 100 • Zb- 


AEW(v.p) 


lim —EgSpE 

^■OO 


2 Ts 


sin 2 {pp 


_(Vp_Y y 

= -Es£pe -sm 2 (^p, 


where is given by 


A snapshot of the simulation box for t = 0 is shown 
in Fig. 5. The use of a Fermi-Dirac distribution instead 
of a Rect-function follows the advice in [46], where it is 
recommended to avoid sharp gradients (which is infinite 
for a Rect-function) because of numerical diffusion and 
spurious oscillations. To ensure the accuracy of the sim¬ 
ulations, we use a sufficiently high number of grid points 
for the Fermi-Dirac-function and generated shock waves 
in order to resolve the gradients properly so that spurious 
effects are suppressed. 

The initial conditions are: 


_ (</>p —^Op)^ 

0Op) ~ 6 COs((^j) 0Op) 

Es (0s, 00s ) — ^s^sR^^Ct ((0s 00s )/^s) 

~ £s^sRD( 0 s 00 s) 

^2 (02 7 0O2 ) - ^2 ^ ^^2(02 7 002) 

with 0i = uJiZ, 0oz = uJiZoi, = uJiTi and i e {p, s}. 

To define the pulse duration Ts of the strong pulse when 
using the Fermi-Dirac function, we equate the calculation 
of the first iteration for the simulational parallel setup 


h^^\(ps) = — f dyFD^{y). 

J-00 


The duration Ts is then defined by Tg = 
lim,^^_^oo The initial conditions are cho¬ 

sen such that the field invariants and the field values at 
the boundary are essentially zero initially and the sys¬ 
tem is simulated until the pulses are again well separated. 


Results of the simulation were compared to the analyt¬ 
ical result for asymptotic lowest order second harmonic 
generation in the parallel and perpendicular set-ups [28], 
where the Gaussian strong background in [28] is replaced 
with the mirrored Fermi-Dirac distribution Eq. (29). 
The excellent agreement is displayed in Fig. 6, where 
the log-log plot of the ratio I{2ujp)/I^\ujp) for various 
values of the strong field amplitude is calculated using: 


4°\ojp) 


\Ep{^)\ 

(^p)l 


Ep(u;) 



(30) 
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from the Schwinger-Dyson equation applied to the 
polarisation operator [52]. 

Photon merging via single four-photon scattering is 
prohibited in a plane wave counterpropagating parallel to 
the background [28, 36]. However, when the possibility 
of multiple four-photon scattering is taken into account, 
high harmonic generation can take place. The modified 
refractive index Eq. (32), experienced by the probe due 
to the strong field and conversely the modified refractive 
index experienced by the strong field due to the probe, 
leads to the electromagnetic invariants Q no longer 
vanishing for the probe and strong fields separately. A 
log-log plot of the normalised spectrum I{uj)/Ip^\ujp) 
for various cases of high harmonic generation through 
purely four-photon scattering is displayed in Fig. 7. 


FIG. 6. (Color online) The relative intensity of the second 
harmonic generated by single six-photon scattering for £p = 
10 "®. 


V. ALL-ORDER FOUR-PHOTON SCATTERING 


For the parameter regime of interest, the most proba¬ 
ble effect on the probe pulse due to four-photon scatter¬ 
ing is that from the well-studied modified vacuum indices 
of refraction = 1-1- given by [50, 51]: 

=2(llT3)/rif2, (31) 


which can be written in a phase-dependent way ni((^s) = 
1 -h 5ni((/Js): 


Sni{(ps)=Am 


4:{ep-esf + 7{epAesf 


(32) 


Following the analytical method in Sec. Ill, summing all 
perturbative orders, one finds due to purely four-photon 
scattering (corresponding to using T = Ti in Eq. (18)), 
in the parallel set-up: 




^pi^Pp) = E 

j=o dPp 


which is just a shift-operator in the phase that is applied 
to the initial probe pulse giving: 

Ep((/3p) = ^ {ipp + ui), 

where the multi-scale parameter for the parallel and per¬ 
pendicular cases t>i = uE: 

up = 2(11 T 3)Mifs $ = (33) 

This all-order solution to the phase shift in a plane wave 
propagating through a constant background derived 
from the Heisenberg-Euler Lagrangian complements 
a recent example solution of the phase shift derived 



uj/cdp 


FIG. 7. (Color online) High harmonic generation from multi¬ 
ple four-photon scattering for ui = 3.3 x 10““*. 


The perpendicular set-up leads to even harmonics be¬ 
ing generated in the mode and odd harmonics being 
generated in the mode. All higher harmonics in the 
perpendicular set-up are suppressed compared to the par¬ 
allel set-up, with odd harmonics being suppressed more 
than even ones. In the parallel set-up, all photons are 
scattered into the mode. As will become clear in Sec. 
VI, compared to the six-photon channel, harmonic gen¬ 
eration via four-photon scattering is considerably sup¬ 
pressed. The scaling argument given at the end of Sec. 
Ill can now be understood in the following way. For 
purely four-photon scattering, one scattering event must 
have occurred to change the electromagnetic variants (a 
factor (5ni = and one further scattering with a 

probe photon (a factor vi = WpiSgSp^), which yields 
the combination 


I/I = (34) 

If one takes this to be the shock parameter for purely 
four-photon scattering, for the parameters of the parallel 
set-up in Fig. 7, following Eq. (20), one would expect the 
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second harmonic at relative intensity (i/i/2)^ = 10“^ ®, 
third harmonic at = 10“^"^ ® and the fourth har¬ 

monic at (z/®/12)^ = 10“^® °, which correctly predict the 
numerical results to within an order of magnitude. For 
comparison, the shock parameter for purely six-photon 
scattering for this set-up would be z /2 = 2.7 z/i. 

VI. ALL-ORDER SIX-PHOTON SCATTERING 

As already hinted, six-photon scattering is the domi¬ 
nant process in the generation of higher harmonics for 
A <g; 1 in the plane wave set-up we are considering. 
For this reason we choose here to analyse six-photon 
scattering as the single vacuum interaction. Many of 
the features of the following harmonic spectra will be 
common to the combined four- and six-photon scattering 
case in Sec. VII. 

The parameter V 2 = ^2 exp[—is bounded by 
V 2 ^ 1 ^ 2 , so the different behaviour of the scattered probe 
will be quantified using the shock parameter U 2 . As V 2 
is increased from zero, two regimes become apparent: i) 
the perturbative regime 1^2 ^ 1 where the occurrence 


of higher harmonics is exponentially suppressed; ii) the 
shock regime, where the intensity of the jth harmonic is 
proportional to a power-law with 7 (v) < —2. 

To highlight the nature of the harmonic generation 
surrounding shock formation, we refer in the following to 
the parallel set-up for simplicity, and discuss differences 
in the perpendicular set-up in Sec. VIC. 

In Fig. 8 are log-log plots of three different types of 
normalised spectrum I{lo) / I^\ujp) in the parallel set-up. 
In the first pane (a) V 2 = 0.05 1 and the perturbative 

regime can be recognised by the exponential suppression 
of higher harmonics. In the middle pane (6) V 2 = 0.6 
and a transition regime can be identified in which the 
lower harmonics are no longer exponentially-suppressed 
but obey a power-law behaviour and the leading-order 
perturbative expansion is inaccurate for higher harmon¬ 
ics. In the final pane (c) V 2 = ^ and the entire plotted 
spectrum has a power-law behaviour, distinctive of the 
shock regime, in which an all-order expansion is required 
to even reach a correct qualitative conclusion. Since we 
are considering only six-photon scattering, we set V 2 = v 
and V 2 = V va. the following discussion. 



FIG. 8. (Color online) Harmonic spectra in the parallel set-up for different regimes of solution: (a) V 2 = 0.05, (fe) U 2 = 
0.6, (c) V 2 = 1. The dots show the leading-order perturbative term, the dashed line is the all-order analytical solution and the 
solid line is from numerical simulation. 


A. Perturbative regime 

li V 1, the amplitude of each harmonic in the scat¬ 
tered electric field is: 


For vj <C 1 but j ^ 1, using Stirling’s approximation 
[53] r(l -h j) Ri we see: 




{y%y 


aj{iy) 


i^e, (36) 


and the exponential dependency of each harmonic 
becomes manifest. In the first pane of Fig. 8, the dots 
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denote the intensities of the harmonics when only the 
leading perturbative order is taken into account. The 
excellent agreement is typical of the perturbative regime, 
in which only a small proportion of probe photons have 
scattered, and double-scattering is much less probable 
than single-scattering. In the transition regime, the 
leading-order terms of the perturbative expansion 
overestimate the intensity of the higher harmonics. In 
the shock regime, the leading-order perturbation terms 
both qualitatively and quantitatively disagree with the 
numerical solution and all-order analytical solution. 


B. Shock regime 



In this regime, v no longer fulfills v ^ 1 and all orders 
of the perturbative expansion must be summed in order 
to calculate the spectrum of generated harmonics. This is 
demonstrated in the third pane of Fig. 8 which shows ex¬ 
cellent agreement between the numerical and analytical 
solution Eqs. (23) and (24). We note that even though 
the all-orders solution includes the phase-dependent pa¬ 
rameter V = z/exp(—((/ 5 p/$p)^), we can still arrive at a 
qualitative understanding of this regime by considering 
the effect on the probe pulse at the point ipp = 0. In this 
case, V = V and the relative amplitude of consecutive 
harmonics is 


FIG. 9. (Color online) After passing through the polarised 
vacuum in the parallel set-up, the probe pulse wavefronts can 
steepen significantly. 

Those parts of the probe that are negative but have a 
larger amplitude are decelerated less than those that are 
negative but have a smaller amplitude, hence leading 
to a steepening in the opposite direction where the 
field is negative. The result is the development of a 
saw-tooth waveform shown in Fig. 9, which is typical of 
a second-order susceptibility [55]. 


aj+i{iy) 


■Ij+ilU + E)v] 

3 

aj{iy) 



j + 1 


(37) 


Using the asymptotic form for jV —>• oo, \Jj{jv)\ 
(27rj)“^/^ (when phase terms are neglected) [54], we see 
that for large enough argument, the ratio of harmonic 
amplitudes becomes: 


\aj{v) \ - 



aj-n(^) 

aj{v) 



(38) 


and the power-law behaviour is manifest. For u = 
1 , this gives a ratio of the intensity of the jth har¬ 
monic to the initial probe intensity, Ip\(pp)/I^\(pp) = 
[Ep\^p) / of 


log 


( Ip'' {^p = 0 ) \ 

\4^\^p = Q)) 


JV—^OO 



31ogj. (39) 


The predicted gradient of 7 = — 3 should be an overes¬ 
timate because for all parts of the probe apart from at 
ipp = 0, V < v. In fact, the full result in the third pane 
of Fig. 8 yields 7 = —3.4. 


A plot of the scattered probe field and induced 
electromagnetic shock is displayed in Fig. 9. Those 
parts of the probe field that are positive and have a 
larger amplitude are decelerated more than those that 
are positive with a smaller amplitude. Where the field 
is positive, this leads to a steepening behind the peaks. 


The coefficient of the jth harmonic is weighted with 
the Bessel function Jj{jv). When v is small, Jj{jv) is 
a rapidly decaying function of j so higher harmonics are 
strongly suppressed. As v —>• 1“, the decay becomes 
much shallower. So a simplified picture of what type 
of shock is generated for the scenario explored in this 
paper can be made by setting the Bessel function to a 
constant. In the parallel set-up, a discontinuous electric 
field with a backwards-leaning waveform of the form Fig. 
9 is generated with: 


E(<^)=ef^(-I)^' 

i=i 


cos( 2 j — \)<p 

27-1 


sin 2 jip 
2 j 


,(40) 


with polarisation e and amplitude £, and the correspond¬ 
ing intensity spectrum has a power law ~ j“^ for har¬ 
monic j. Indeed we find on a plot of 7 )^) (see Fig. 10), 
that as v increases above 1 , the power-law exponent in 
the numerical spectrum increases, tending toward a the¬ 
oretical maximum of — 2 , at which point the lack of a 
unique solution to Maxwell’s equations would halt fur¬ 
ther propagation of the probe. For v > 1, the numerical 
spectrum displays a variable power law, which is shal¬ 
lower for higher harmonics where the agreement with the 
analytical solution Eqs. (23) and (24) becomes increas¬ 
ingly worse. The power law exponent calculated using 
the fourth and tenth harmonic is displayed in Fig. 10, 
where unlike in the numerical solution, in which the spec¬ 
trum becomes progressively shallower, the analytical so¬ 
lution reaches a maximum shallowness. It is unclear what 
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physical mechanism would cause this maximum to occur, 
which suggests this is a limitation of the viability of the 
analytical solution. Indeed when u > 1 in the analytical 
solution, Jj{jv) can oscillate with j, and the ordering of 
harmonics can become no longer monotonic. 


As the numerical spectrum becomes shallower, very high 
harmonics appear, which questions the validity condition 
jujp <C m for using the Heisenberg-Euler Lagrangian to 
describe the vacuum interaction, and questions how steep 
the power law can become before relaxation processes 
would take over. 


- 2.0 

-2.5 

-3.0 

-3.5 

-4.0 

-4.5 

0.8 0.9 1.0 1.1 1.2 1.3 1.4 

v 

FIG. 10. (Color online) Comparison of the power-law expo¬ 
nent ^iy) for different values of the shock parameter v, as 
calculated using the fourth and tenth harmonics from the an¬ 
alytical (dashed line) and numerical (points) solutions. 



C. Polarisation Dependency 

The previous sections are for the parallel set-up. For 
the perpendicular set-up, even harmonics are generated 
in the parallel mode and odd harmonics in the 

perpendicular mode This is demonstrated in the 

spectrum in Fig. 11, where the thick and thin lines 
distinguish how the generated harmonics are polarised. 

The shock wave generated in the perpendicular set¬ 
up is displayed in Fig. 12. The scattered field in the 
mode demonstrates a shock of a different nature to in the 
parallel set-up, tending towards a square rather than a 
saw-tooth waveform. Such a waveform can be generated 
with the sum: 

= <=£ £(- 1 )’ (41) 

1=1 ^ 


which is just the odd frequencies of Eq. (40). 
J 



3 



FIG. 11. (Color online) Harmonic spectra from numerical simulation of the perpendicular set-up for different regimes of solution: 
(a) V 2 = 0.05, (6) V 2 = 0.6, (c) V 2 = 1. The thick blue (thin green) peaks are harmonics parallel to the probe (strong) pulse. 


In the mode, a similar shock to in the parallel set-up 
is seen, only with double the frequency. Such a saw-tooth 


electric field is given by the sum [53]: 


1=1 


sin 2 jip 


( 42 ) 
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which is just the even frequencies of Eq. (40), beginning 
at double the frequency of the seed probe field. 



FIG. 12. (Color online) A probe that is initially polarised per¬ 
pendicular to the background (blue dashed line) experiences 
different shocks in the (dot-dashed red line) and (solid 
green line) modes (color online). 


VII. ALL-ORDER FOUR- AND SIX-PHOTON 
SCATTERING 

Although six-photon scattering is the most efficient 
process in generating high harmonics, for the parameter 


regime we are interested in, the effect of four-photon 
scattering as a modified vacuum refractive index cannot 
be neglected. Since the interaction with the vacuum 
includes powers of the probe field, the effects of phase lag 
and harmonic generation can mix in a highly nonlinear 
way. In this section we give the results of numerical 
simulations that include both processes. 


For the parallel set-up, the spectrum generated by six- 
photon scattering (for example, as shown in Fig. 8), is 
not visibly affected by the inclusion of four-photon scat¬ 
tering. However, for the perpendicular set-up, since even 
and odd harmonics are in different polarisation modes 
and since the vacuum is birefringent so each polarisa¬ 
tion mode experiences a different phase lag, the inclu¬ 
sion of four-photon scattering was found to increase the 
asymmetry between the even and odd harmonics com¬ 
pared with the purely six-photon scattering case. This is 
demonstrated in Fig. 13 for the case vi = 100, 1/2 = 1, 
which compares the spectrum of harmonics generated 
when: i) only four-photon scattering is included (left- 
hand pane); ii) only six-photon scattering is included 
(middle pane) and hi) four- and six- photon scattering are 
included (right-hand pane). The right-hand pane demon¬ 
strates the increased asymmetry between even and odd 
harmonics. 



FIG. 13. (Color online) High harmonic generation for the perpendicular set-up when four- and six-photon scattering are present 
and four-photon scattering is much more prevalent than six-photon scattering (wi = 100, V 2 — 1). The first pane (a) is for just 
four-photon scattering, the second pane (6) for just six-photon scattering and the third pane (c) for when both are present. 
The thick blue (thin green) peaks are again harmonics parallel to the probe (strong) pulse. 


As the case of four- and six-photon scattering differs from the six-photon scattering case only for the perpen- 
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dicular set-up, we focus our discussion on this. Then 
there are three cases of interest: i) weak dispersive: 
vi <C 1 ^ 2 ', ii) dispersive: vi ~ V 2 ] ih) strong dispersive: 
v\ ^ V 2 . The first case of weak vacuum dispersion is 
within the parameter regime of interest, but outside of 
the regime that can be numerically simulated as it would 
require ^ 12 ^ > 0.1 if the hierarchy £s » £p were to be 
maintained. In the limit of vanishing dispersion, we ex¬ 
pect the results from purely six-photon scattering case to 
be valid (this will be seen to be implied from the results 
of a dispersive vacuum). 

A. Dispersive vacuum vi ~ V2 

When vacuum dispersion is significant, one might 
expect the nature of the shock wave to change. Two 
cases were simulated: i) when ui = z /2 = 1 and ii) when 
Ui = 5, z /2 = 1. For the first case of equal parameters, 
the shock wave in Fig. 14 was generated. This bears 
a close resemblance to the shock wave generated in 
the perpendicular set-up for a dispersionless vacuum 
{v\ 0), i.e. when only six-photon scattering is present, 

but with a noticeable lag due to the now non-unitary 
refractive index. However, when the amount of dis- 



FIG. 14. (Color online) The weakly-dispersive case for 
the perpendicular set-up. The initially e'l polarised probe 
(blue dashed line) experiences the mixture of the probe- 
independent vacuum refractive index (here t>i = 1) and 
the shock-inducing probe-dependent vacuum refractive index 
(here V 2 = 1). The mode (dot-dashed red line) and 
mode (solid green line) behave differently. 

persion is increased, setting ui = 5 and V 2 = 1, the 
shock wave takes on the different form shown in Fig. 
15. In this dispersive case, the parallel mode develops 
a shock reminiscent of an optical Kerr medium, in 
which the polarisation contains a cubic nonlinearity 
Pi = + xifkiPjPkEi. This is in some ways 

unsurprising because the parallel mode only contains 
odd harmonics and therefore odd powers of the held. 



FIG. 15. (Color online) When the dispersion is increased 
(ui = 5, 1^2 = 1), a different type of shock is formed in the 

mode (dot-dashed red line) and the shock in the mode 
(solid green line) is reduced, where the initially e'' polarised 
probe is plotted by the blue dashed line. 

and the largest nonlinear term originates from an Ep 
term. Therefore the symmetry of the scattered held 
when the held direction is swapped £p —>■ —£p is different 
for the parallel held (which contains only even powers 
of £p) and the perpendicular held (which contains only 
odd powers oi £p). 

Carrier-wave shocking also occurs in nonlinear opti¬ 
cal materials. Our hndings are similar to those reported 
in [56], where excellent agreement was obtained between 
theory and simulation in the dispersionless limit of a 
Kerr-like nonlinear material, but where it was noted how 
involved the analysis becomes if there is a complicated 
phase dependency between the generated harmonics. In 
the current work, in the parallel set-up with dispersion 
(i.e. four- and six-photon scattering present), all harmon¬ 
ics experience the same refractive index so a shock wave 
can build up. In the perpendicular set-up, the refractive 
index in the mode is different to in the mode. We 
are studying a regime in which harmonics are generated 
by a chain of scattering processes. Since, in each chain 
of processes that lead to the generation of a specific har¬ 
monic, the probe spends a different amount of time in the 
than in the mode, the probability for each chain 
will be multiplied by a different phase. When the proba¬ 
bility of all possible chains is summed over, it is reduced 
compared to the parallel set-up due to each probability 
being added incoherently. This leads to a suppression of 
shock wave generation. 


B. Strongly-dispersive vacuum vi ':$> V2 

To investigate shock wave generation in the strongly- 
dispersive regime, we set vi = 100 and j ^2 = 1- A new 
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type of behaviour becomes apparent, namely the defor¬ 
mation of the probe pulse envelope. The bandwidth of 
the probe is of the order 1 /r^ but due to dispersive effects, 
frequencies of this magnitude can no longer be neglected. 
Since vi = Sipp = tOpT, where T is the duration of prop¬ 
agation, frequencies from the probe envelope separated 
by 1 /Tp will acquire a temporal separation relative to the 
duration of the pulse of vijuipTp 1. Furthermore, the 
second harmonic is considerably suppressed when disper¬ 
sion is included, such that it is of the same order of mag¬ 
nitude as the scattering of the probe envelope frequency. 
For this reason, the effect on the probe envelope can be 
seen so clearly in the £"*■ component in Fig. 16. 



FIG. 16. (Color online) The probe pulse after having scat¬ 
tered in the strong background when ui = 100, 1^2 = 1. 


VIII. DISCUSSION 

A. Comparison with high harmonic generation in 
oscillating plasmas 

There is a certain similarity between high harmonic 
generation due to the relativistic movement of electrons 
in the plasma of laser-irradiated foil experiments and the 
virtual electron-positron “plasma” of the laser-irradiated 
vacuum. The vacuum is transparent when the invariants 
£2 - ^2 and E • B are zero. Therefore the vacuum is 
transparent to a pure plane wave and these invariants 
also typically remain much smaller for a single focused 
pulse than for counterpropagating pulses. So unlike with 
the plasma present in a foil, the vacuum plasma is first 
“activated” by being polarised by some second “pump” 
pulse, similar to in a pump-probe experiment. In the 
current work, the vacuum was polarised by a back¬ 
ground with the profile of a rectangular function. As the 
leading-order nonlinear polarisation was proportional to 
the applied field cubed, it suggests that the local charge 
density is also non-zero in this region. The rectangular 
function is used to model the electron density in a solid 


before it is exposed to a laser pulse [57] and also to 
represent the laser’s profile and in capillary discharge 
waveguides [58]. The difference with the vacuum is that 
the polarised material can in some way be “formed” 
by the pump pulse in the moment it is traversed by a 
probe. 

For the parallel set-up, all harmonics were generated 
in the parallel mode, but for the perpendicular set-up 
odd harmonics were generated in the parallel mode 
with even harmonics in the perpendicular one. Just 
as in single nonlinear Compton scattering [59], the 
generation of the parallel mode is more probable than 
the perpendicular one. The relationship between polar¬ 
isation and harmonic order is reminiscent of selection 
rules for harmonics generated in laser-foil experiments, 
for example in the “p-polarised” (parallel to plane of 
incidence) and “s-polarised” (perpendicular to plane 
of incidence) harmonics in the widely-used oscillating 
mirror model [60]. 

In the harmonic spectrum generated by a real plasma 
in laser-foil experiments, there is also a region of 
power-law decay and a region of exponential decay, 
as found here for vacuum high harmonic generation. 
For the oscillating mirror model, power-law exponents 
of 7 = -5/2 [61] and 7 = - 8/3 [ 52 ] have been postu¬ 
lated, and experiments on solid targets have recorded 
intensity-dependent power-law exponents, for example 
in [63] of —5.50 < 7 < —3.38. These values are close to 
the analytical and numerical values found in the current 
work for vacuum high harmonic generation in the shock 
regime, —4.5 ^ 7 < —2 (the lower limit corresponds to 
the gradient when the power-law behaviour becomes 
manifest at 1^2 ~ 0.85). Moreover, the power-law expo¬ 
nent 7 = ^{ 1 ^ 2 ) is also a function of the shock parameter 
1^2 oc $ and therefore increases with further propagation 
of the probe through the polarised vacuum, up to a 
theoretical maximum of 7 ( 1 ^ 2 ) < —2. In contrast to 
the overdense plasma case, with our plane-wave model 
and increasing shock parameter, we found no indication 
of a frequency cutoff, although at some frequency, 
pair-creation processes will play a role. By this we mean 
that higher harmonics can seed tunnelling pair creation 
in the background field [64, 65] or colliding photons with 
wavevectors fci and k 2 satisfying kik 2 > 2 to^ would lead 
to multi-photon pair creation (the Breit-Wheeler process 
[ 66 , 67]). This would presumably deplete the higher 
harmonics that are directly related to the steepening of 
the wave fronts and act as a wave-breaking mechanism 
for the shock wave. 

In the parallel set-up, each harmonic has a regular 
phase relationship to the others and so a shock wave can 
build up as the amplitude of higher harmonics increases. 
In contrast to this, in the perpendicular set-up, since 
there are many different chains of processes that can 
lead to the creation of a given harmonic, and since in 
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each chain a different amount of time is spent in each 
polarisation, which leads to different dispersion relations, 
the phase of each harmonic is related to the others in a 
non-trivial way and they are summed incoherently. This 
behaviour is similar to that found in studies of non-linear 
optical materials [56], and leads to the suppression of 
shock wave generation. 

Although high harmonic generation is present in laser- 
gas and laser-liquid experiments, the spectrum generated 
is of a completely different form. As harmonics are gen¬ 
erated via the three-step recombination mechanism un¬ 
dergone by an electron in the Coulomb field of a nucleus, 
the electrons’ trajectory and hence harmonics generated, 
are of a fundamentally different nature and demonstrate 
a genuine “plateau” region in spectra that is not present 
in vacuum high harmonic generation as studied in the 
present work [ 68 ]. 


B. Validity of approach 

By considering colliding plane waves, scattering in the 
transverse direction was ignored. One can estimate when 
this is a good approximation by defining the diffraction 
parameter I = w'^/XpTs, where w is the width of the 
probe pulse in the transverse plane (assumed smaller 
than the width of the background). When I ^ 1 one 
is in the “near zone” and diffraction effects should be 
negligible whereas I <C 1 represents the “far zone” and 
diffraction effects become important [69]. 

The numerical simulation and analytical calculation 
predict a self-steeping of the probe wavefronts, which 
increases with shock parameter iz, until the wavefronts 
reach a theoretical maximum of becoming infinitely steep 
at which point the solutions to the wave equation are no 
longer unique. Since the Heisenberg-Euler Lagrangian 
is expected to be valid when the typical scale of a field 
inhomogeneity is much larger than the reduced Compton 
wavelength, this infinite steepening is not expected to be 
physically realisable. Moreover, no relaxation processes 
are included. If transverse dimensions would be taken 
into account, since six-photon scattering depends on 
the probe amplitude, self-focusing effects should be 
present. Furthermore, self-focusing can also occur via 
four-photon scattering as the probability for asymptotic 
second-harmonic generation via four-photon scattering 
becomes non-zero when the colliding probe photons do 
not propagate in parallel. So when transverse dimensions 
are included, as the probe propagates, it becomes less 
like a plane wave and the higher harmonics can seed real 
electron-positron pair creation as previously described. 

The polarisation of other vacuum virtual particle 
species such as muons, pions and quarks was neglected, 
as the energy scale associated with these particles is much 
higher [70]. For that reason, we confined our discussion 


to the polarisation of virtual electron-positron pairs. 


C. Measurability 

Vacuum high harmonic generation in the shock 
regime becomes important when the shock parameter 
zz Ri 1. Taking as an example six-photon scattering 
for the parallel set-up, u = V 2 = The 

current record for the highest electric field of a laser 
pulse produced in a laboratory [71] is of the order 
3 X 10“^Acr. Recalling that fields are written in units 
of the critical field, and that ^2 = <C 1 , it is 

clear that the shock regime is currently well out of 
the reach of optical laser-based experiments. Vacuum 
polarisation effects that can more likely be measured in 
laser-based experiments include elastic photon-photon 
scattering [29, 72-83] or lowest-order photon merging 
[26, 27, 84, 85]. (A review of strong-held QFD effects 
can be found in [ 86 , 87].) The current best experimental 
limits for photon-photon scattering in an all-optical laser 
set-up [ 88 ] and combining magnetic helds with resonant 
optical cavities [89, 90] are still orders of magnitude 
above the QFD prediction. 

Where such vacuum electromagnetic shocks and ac¬ 
companying harmonic generation might play a role, is in 
the evolution of X-ray pulsars and strongly-magnetised 
neutron stars or “magnetars” [91-93]. Photons are emit¬ 
ted from the surface of such objects and propagate 
through magnetic helds of strength up to and beyond 
(in the system of units we use, Rcr = I^cr), in plasmas of 
around 0.1-10cm in depth [94]. The current results were 
derived for a constant crossed held background, but can 
be generalised to a constant magnetic held, which should 
be a good approximation to the local held in strongly- 
magnetic pulsars, which is expected to be that of a dipole 
[94] on the stellar scale. 

IX. CONCLUSION 

When the quantum nature of the vacuum is taken 
into account, an electromagnetic shock accompanies high 
harmonic generation in an oscillating plane probe pulse 
counterpropagating through a stronger slowly-varying 
plane pulse. We have identihed a nonlinear shock pa¬ 
rameter that indicates when the self-interaction of the 
probe due to the polarised vacuum becomes important 
and have shown that this can be consistently described 
using a probe-dependent vacuum refractive index. 

As the shock parameter increases from zero, the spec¬ 
trum of generated harmonics moves from an exponential 
decay to a power-law decay. The intensity of the jth har¬ 
monic in the shock regime was found in an all-order ana¬ 
lytical solution and numerical simulation to be p, where 
7 increases with propagation distance. A power law be¬ 
haviour was observed for —4.5 ^ J ~2.4, where the 
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exponent is theoretically limited by 7 < —2 as the probe 
pulse wavefronts would become infinitely steep and could 
no longer propagate. Due to the very high generated 
frequencies, the Heisenberg-Euler approach is no longer 
applicable at this point. Moreover, relaxation processes 
such as photon-seeded and Breit-Wheeler pair creation 
should then become probable. 

When the polarisation of the probe and background 
is parallel, all harmonics are generated in the parallel 
mode, but when the probe is perpendicularly-polarised 
to the background, odd and even harmonics are split 
into probe and background polarisation modes respec¬ 
tively. Due to the birefringence of the vacuum, the probe 
polarisation mode is generated more abundantly than in 
the background polarisation mode. Moreover, due to the 
separation of frequencies, the parallel set-up displays a 
saw-tooth shock in the parallel mode, whereas the per¬ 


pendicular set-up displays a Kerr-like shock. 

Both the simulational and analytical methods pre¬ 
sented can be generalised to more complicated probe and 
background helds. 
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Appendix A: Coefficients for modified Maxwell 
Equations 


We 

define Cxy = d^LuE/dxdy and 

r^ = a? -\- b^. 


Cl 

dHa bILi) 
= 47r- — - 

5 



(Al) 

C 2 

= TT^ [a(a^ - 36^)£a 

+ - 

Za^)Cb 



{a^Caa 

, — 2 ab Cab + b^ 

Cbb)] 

(A2) 

C 3 

= - 

a^)Ca 

-1- 6(3a^ - 

- b^)Cb 



{b^Caa 

+ 2 a 6 Hah 

Cbb)] 

(A3) 

C 4 

= - 


-1- a{a^ — 

ib'^)Cb 



(abCaa 

+ (a 2 

- b^) Cal 

■} dh Hh}^^ 

(A4) 


For the first order (box diagram) and the second order 
(hexagon diagram) in the weak-field expansion, we find 
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the following coefficients: 


Appendix C: Biased finite differences 


C^l,Box 

C^3,Box 


457r^ ^ 

^C'2,Box 


C 2 ,Box = - 7 ^ (A5) 

457r 

C'4.Box = 0 (A 6 ) 


Ci.Hex = - B^f + 13(E • Bf] (A7) 

ol^TT 

A.f\i 

C 23 ex = -^(A^-S^) (AS) 

C3,Hex = ^(1^2,Hex (A9) 

1 Srv 

C4,Hex = - 7 T^|(E-B)| (AlO) 

3157r 


The action of the matrix D on the vector u can be 
encoded in the use of an adaption of the DSS020 function 
from [46]: 


Du = 


/d-{u\)\ 

d-lul) 

d+{ul) 

d+{ul) 

d-luf) 

d-lul) 

d+{ul) 

d+{ul) 


V ; / 


Appendix B: Matrix form of modified Maxwell 
Equations 

The modified Maxwell equations Eqs. (5) and ( 6 ) can 
be written in matrix form: 

(l 4 + X)atf+(Q + Y)5,f = 0, (Bl) 

where f = [Ex, Ey, Bx, By)"’", I 4 is the identity in four 
dimensions, Q = adiag(l, —1, —1,1) and X = (xij), Y = 
(Uij) are the vacuum perturbation, where the non-zero 
elements are given by: 

2^11 = Cl — C2P11 — C3P33 — 2C4P13 
a;i2 = —C2P12 — C3P34 — C4{pi4 -I- P23) 

2:13 = iC2 — C3)pi3 + C4{p33 — Pii) 

Xi 4 = C 2 P 14 — C 3 P 23 + C'4(p34 — P 12 ) 

2:21 = —C 2 P 12 — C 3 P 34 — C 4 {pi 4 -f P 23 ) 

X 22 = Cl — C 2 P 22 — C 3 P 44 — 2 C 4 P 24 
2:23 = C2P23 — C3P14 + C 4 {p 34 — P12) 

X 24 = {C 2 — C3)p24 + C 4 {p 44 — P 22 ) 
yil = —C 2 P 14 + C 3 P 23 + C4(pi2 — P 34 ) 
yi2 = —{C2 — C 3 )p 24 + C 4 (p 22 — P 44 ) 
yi3 = C2P34 + C3P12 — 64(^14 + P23) 
yi 4 = Ci+ C2P/14 + C3P22 — 2C4P24 

?/21 = {C 2 — C3)pi3 + C4{p33 — pii) 

2/22 = C2P23 — C3P14 -f C4(p34 — P12) 

2/23 = —Cl — C 2 P 33 — C 3 P 11 -f 2 C 4 P 13 
2/24 = —C 2 P 34 — C 3 P 12 -f C 4 (pi 4 -f P 23 ) 

where we define pij := 4:fifj, such that e.g. pi 4 = AExBy. 


where the function d-{u’‘) is defined as 

d-{u’) := 

1 = 1 : 

q{—25u’ + 48 m^ — 36 u^ -I- — 3 m^) 

l = N -2: 

-b 8u^-i - 

l = N -1: 

-b - 18 m^“ 2 -b -b 3 m^) 

l = N: 

q{3u’^-'’ - -b 36 u^-2 _ + 25m^) 

else : 

q{-3u^~’ - 10m' -b 18 m'+^ - 6m'+^ -b 

with q = 1 / 12 Az and d+{u’) as 

d+(M') := 

1 = 1: 

q{—2hu’ + 48u^ — + 16u^ — 3m^) 

1 = 2: 

q{—3u’ — lOu^ -b ISm"^ — 6m"' -b M®) 

1 = 3: 

q{u’ - 8m^ -b 8 m"' - m") 

l = N: 

q{3u’^-'’ - + 36m'^- 2 _ 48uAf-i + 25 m'^) 

else : 

qi-u’-^ + 6m'-2 _ 18m'-' -b 10 m' -b 3 m'+'). 







